#load gplot2 library
library(ggplot2)
#Read data 
library(readxl)
Data <- read_excel("~/Charles/Research/Phytoremediation of Fe3+ soil/Data/Data.xlsx", 
                    sheet = "Iron in soil")
View(Data)
#Remove rows with NAs
Data <- na.omit(Data)
View(Data)
#convert categorical variables to factors 
str(Data)
Data$Soil_mixture = as.factor(Data$Soil_mixture)
Data$DAT = as.factor(Data$DAT)
Data$Planted_Unplanted = as.factor(Data$Planted_Unplanted)
Data$Plant_height = as.numeric (Data$Plant_height)
Data$Leaf_length = as.numeric (Data$Leaf_length)
Data$Leaf_width = as.numeric (Data$Leaf_width)
Data$Iron_soil = as.numeric (Data$Iron_soil)
Data$Percent_increase_plant_height = as.numeric (Data$Percent_increase_plant_height)
Data$Percent_increase_leaf_length = as.numeric (Data$Percent_increase_leaf_length)
Data$Percent_increase_leaf_width = as.numeric (Data$Percent_increase_leaf_width)
Data$Leaf_biomass = as.numeric (Data$Leaf_biomass)
Data$Stem_biomass = as.numeric (Data$Stem_biomass)
Data$Leaf_biomass = as.numeric (Data$Leaf_biomass)
Data$Root_biomass = as.numeric (Data$Root_biomass)
Data$Total_biomass = as.numeric (Data$Total_biomass)
Data$Root_shoot_ratio = as.numeric (Data$Root_shoot_ratio)
Data$Iron_leaf = as.numeric (Data$Iron_leaf)
Data$Iron_stem = as.numeric(Data$Iron_stem)
Data$Iron_roots = as.numeric(Data$Iron_roots)
Data$BCF_leaf = as.numeric(Data$BCF_leaf)
Data$BCF_stem = as.numeric(Data$BCF_stem)
Data$BCF_root = as.numeric(Data$BCF_root)
Data$TF_leaf = as.numeric(Data$TF_leaf)
Data$TF_stem = as.numeric(Data$TF_stem)
str(Data)
View(Data)


#two-way ANOVA with no interaction 
two.way <- aov(Iron_soil ~ Soil_mixture + Planted_Unplanted + DAT, data = Data)
#two-way ANOVA with interaction
interaction <- aov(Iron_soil ~ Soil_mixture * Planted_Unplanted * DAT, data = Data)
#Levene's test
install.packages('car')

library(car)
leveneTest(TF_stem ~ Soil_mixture, data = Data)

#one-way ANOVA 
one.way <- aov(TF_stem ~ Soil_mixture, data = Data)
#Summary of one-way ANOVA
summary(one.way)
#AIC calculates the best-fit model by finding the model that explains the 
#largest amount of variation in the response variable while using the fewest 
#parameters (Bevans, 2023). 
#https://www.scribbr.com/statistics/two-way-anova/#:~:text=A%20two%2Dway%20ANOVA%20is,combination%2C%20affect%20a%20dependent%20variable.
#install the AIC package
install.packages("AICcmodavg")
#Load the AIC package
library(AICcmodavg)
#Run the AIC package
model.set <- list(two.way, interaction)
model.names <- c("two.way", "interaction")

aictab(model.set, modnames = model.names)


#Summary of two-way/interaction ANOVA
summary(interaction)


#Check for homoscedasticity. 
#Check whether the model fits the assumption of homoscedasticity
par(mfrow=c(2,2))
plot(one.way)
par(mfrow=c(1,1))

#Post-hoc test. 
#Tukey’s Honestly-Significant-Difference (TukeyHSD) test
turkey.two.way<-TukeyHSD(two.way)
turkey.two.way
#Tukey’s Honestly-Significant-Difference (TukeyHSD) test
turkey.interaction<-TukeyHSD(interaction)
turkey.interaction
options (max.print = 1000000)
#Tukey’s Honestly-Significant-Difference (TukeyHSD) test
turkey.one.way<-TukeyHSD(one.way)
turkey.one.way
#plot of turkey test
tukey.plot.aov<-aov(Leaf_width ~ EMS:Accession:DAP, data = Data1)
tukey.plot.test<-TukeyHSD(tukey.plot.aov)
tukey.plot.test
plot(tukey.plot.test, las = 1)

#Presentation of result in bar chart
TF_stem <- ggplot(Data, aes(x = Soil_mixture, y = TF_stem    
                                                     
                                                    
                                                   
) 
) +
  
  
    stat_summary(fun = "mean", geom = "bar", 
               position = "dodge", 
               colour = "black", 
               alpha = .7)+
  stat_summary(fun = "mean", 
               geom = "point",
               position = position_dodge(0.9),
               size = 1)+
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  
  stat_summary(fun.data = 'mean_se', geom = 'errorbar', na.rm = TRUE,  
               position = position_dodge(0.9), width = 0.4, color = 'black')+
  theme(
    strip.text.x = element_text(margin = margin(2, 0, 2, 0))
  )+
  ylim(0, 2) +
  scale_fill_grey()
TF_stem    
TF_stem   <- TF_stem   +
  labs(title = "Translocation factor of iron in stem in response to rusty soil mixture",
       x = "Soil mixture",
       y = "Bioconcentration factor") +
  labs(fill = "Days 
after
transplant
(DAT)
 ")+
  theme(axis.title = element_text(color = "black",
                                  face='bold',
                                  size = 12), 
        legend.title=element_text(color = "black", 
                                  face = "bold", 
                                  size = 12),
        legend.text = element_text(color = "black", 
                                   face = "bold", 
                                   size = 12), 
        plot.title = element_text(color = "black", 
                                  face = "bold", 
                                  size = 12), 
        axis.text = element_text(color = "black",
                                 face = "bold",
                                 size = 8)) +
  geom_hline (yintercept = 1, color = "black", linetype = "solid")

TF_stem   
ggsave("TF_stem.png")

